#include <stdio.h>
#include <math.h>

/**
 *     a function for computing bivariate t probabilities.
 *
 * mvbvtl - calculate the probability that x < dh and y < dk.
 *
 * parameters
 *
 *   nu number of degrees of freedom
 *   dh 1st lower integration limit
 *   dk 2nd lower integration limit
 *   r   correlation coefficient
 *
 */
double mvbvtl (int nu, double dh, double dk, double r) {

	int j, hs, ks;
	double tpi, pi, ors, hrk, krh, bvt, snu;
	double gmph, gmpk, xnkh, xnhk, qhrk, hkn, hpk, hkrn;
	double btnckh, btnchk, btpdkh, btpdhk, one;

	//parameter ( pi = 3.14159265358979323844d0, tpi = 2*pi, one = 1 )
	// Initialisation des valeurs
	tpi = 2*pi;
	one = 1;
	snu = sqrt( double(nu) );
	ors = 1 - r*r;
	hrk = dh - r*dk;
	krh = dk - r*dh;

	if ( abs(hrk) + ors > 0 ) {
		xnhk = (hrk^2)/( (hrk^2) + ors*( nu + (dk^2) ) );
		xnkh = (krh^2)/( krh^2) + ors*( nu + (dh^2) ) );
	}
	else{
		xnhk = 0;
		xnkh = 0;
	}

	hs = sign( one, dh - r*dk );
	ks = sign( one, dk - r*dh );

	if ( nu%2 == 0 ) {
	         bvt = atan2( sqrt(ors), -r )/tpi;
	         gmph = dh/sqrt( 16*( nu + (dh^2) ) );
	         gmpk = dk/sqrt( 16*( nu + (dk^2) ) );
	         btnckh = 2*atan2( sqrt( xnkh ), sqrt( 1 - xnkh ) )/pi ;
	         btpdkh = 2*sqrt( xnkh*( 1 - xnkh ) )/pi ;
	         btnchk = 2*atan2( sqrt( xnhk ), sqrt( 1 - xnhk ) )/pi  ;
	         btpdhk = 2*sqrt( xnhk*( 1 - xnhk ) )/pi ;
	         for (j=1; j<nu/2; j++) {
	            bvt = bvt + gmph*( 1 + ks*btnckh ) ;
	            bvt = bvt + gmpk*( 1 + hs*btnchk ) ;
	            btnckh = btnckh + btpdkh  ;
	            btpdkh = 2*j*btpdkh*( 1 - xnkh )/( 2*j + 1 )  ;
	            btnchk = btnchk + btpdhk;
	            btpdhk = 2*j*btpdhk*( 1 - xnhk )/( 2*j + 1 );
	            gmph = gmph*( 2*j - 1 )/( 2*j*( 1 + (dh^2)/nu ) );
	            gmpk = gmpk*( 2*j - 1 )/( 2*j*( 1 + (dk^2)/nu ) );
	         }
	}
	      else {
	         qhrk = sqrt( (dh^2) + (dk^2) - 2*r*dh*dk + nu*ors ) ;
	         hkrn = dh*dk + r*nu  ;
	         hkn = dh*dk - nu  ;
	         hpk = dh + dk ;
	         bvt = atan2(-snu*(hkn*qhrk+hpk*hkrn),hkn*hkrn-nu*hpk*qhrk)/tpi;
	         if ( bvt < -1E-15 )
	        	 bvt = bvt + 1;
	         gmph = dh/( tpi*snu*( 1 + (dh^2)/nu ) )  ;
	         gmpk = dk/( tpi*snu*( 1 + (dk^2)/nu ) )  ;
	         btnckh = sqrt( xnkh )  ;
	         btpdkh = btnckh ;
	         btnchk = sqrt( xnhk )  ;
	         btpdhk = btnchk  ;
	         for (j=1; j < ( nu - 1 )/2; j++) {
	            bvt = bvt + gmph*( 1 + ks*btnckh ) ;
	            bvt = bvt + gmpk*( 1 + hs*btnchk ) ;
	            btpdkh = ( 2*j - 1 )*btpdkh*( 1 - xnkh )/( 2*j );
	            btnckh = btnckh + btpdkh  ;
	            btpdhk = ( 2*j - 1 )*btpdhk*( 1 - xnhk )/( 2*j ) ;
	            btnchk = btnchk + btpdhk  ;
	            gmph = 2*j*gmph/( ( 2*j + 1 )*( 1 + (dh^2)/nu ) );
	            gmpk = 2*j*gmpk/( ( 2*j + 1 )*( 1 + (dk^2)/nu ) );
	         }
	      }

}



// Renvoie le signe de a*b
int sign(int a, int b) {

	if ( a*b > 0)
		return 1;
	else return -1;

}
